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Abstract. We have characterized the physical structure and chemical composition of two close-to-round starless 
cores in Taurus- Auriga, L1498 and L1517B. Our analysis is based on high angular resolution observations in at 
least two transitions of NH3, N2H+, CS, C'^'*S, C'*0, and C'^0, together with maps of the 1.2 mm continuum. 
For both cores, we derive radial profiles of constant temperature and constant turbulence, together with density 
distributions close to those of non-singular isothermal spheres. Using these physical conditions and a Monte Carlo 
radiative transfer model, we derive abundance profiles for all species and model the strong chemical differentiation 
of the core interiors. According to our models, the NH3 abundance increases toward the core centers by a factor of 
several (~ 5) while N2H"'" has a constant abundance over most of the cores. In contrast, both C'*0 and CS (and 
isotopomers) are strongly depleted in the core interiors, most likely due to their freeze out onto grains at densities 
of a few lO** cm~'^. Concerning the kinematics of the dense gas, we find (in addition to constant turbulence) a 
pattern of internal motions at the level of 0.1 km s~^. These motions seem correlated with asymmetries in the 
pattern of molecular depletion, and we interpret them as residuals of core contraction. Their distribution and size 
suggest that core formation occurs in a rather irregular manner and with a time scale of a Myr. A comparison of 
our derived core properties with those predicted by supersonic turbulence models of core formation shows that 
our Taurus cores are much more quiescent than representative predictions from these models. In two appendices 
at the end of the paper we present a simple and accurate approximation to the density profile of an isothermal 
(Bonnor-Ebert) sphere, and a Monte Carlo-calibrated method to derive gas kinetic temperatures using NH3 data. 
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1. Introduction 

Dense cores like those in Taurus and other nearby dark 
clouds represent the simplest environments where stars 
form. Their study holds the promise of revealing the basic 
physics of star formation, but despite their apparent 
simplicity, dense cores remain mysterious in their internal 
structure, formation and evolution, and in t he manner in 
which t hey collapse under gravity (see, e.g.. iMverslllQQgt 
lEvans I Q jj99 for recent reviews) . Part of this mystery 
results from a number of contradictions betw een observa- 
tions of different dense gas tracers (see, e.g.. IZhou et al. I 
Il989l) . which have so far precluded a self consistent 
description of the core internal properties. Fortunately, 
multi-molecular and continuum observations over the last 
few years have started to suggest that most contradictions 
arise from cores not being chemically homogeneous, as ini- 
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Understanding and characterizing the chemical structure 
of dense cores has therefore become a top priority for star 
formation studies, as solving this problem may finally 
produce a coherent picture of dense cores, from which 
their star forming evolution may be inferred. 

If cores have strong chemical composition gradients, in 
addition to density gradients, their description cannot rely 
on average parameters, as often assumed. It is necessary 
to use spatially- variable abundances, and these need to be 
derived from the data through a realistic modeling of the 
radiative transfer. This approach, unfortunately, involves 
a number of uncertainties, one of the largest being our 
ignorance of the exact (3D) geometry of a core, which is 
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not well constrained either by the ory or 2 D images (e.g., 
compare IJones fc Basu1l2002l with lCurry| |2002). To min- 
imize this effect, it is convenient to concentrate on cores 
with relatively simple geometries, and for this reason, we 
have selected L1498 and L1517B as subjects of the present 
study. 

LI 498 and L1517B are two dense cores in the 
Ta urus- Aurig a cloud complex, and were first identified 
by iMvers et al. I l|l98.1) from the inspection of Palomar 
Sky Atlas prin ts. They are not associated with IRAS 

or 1.2 millime- 
and can there- 
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differentiated system, and furth er details of its chemica l 
structure ha ve been presented bvlWolkovitch et al. 1 ljl997ll 
for CCS andlWil1a,cy et aL I (|l99S^ for CO. More recently. 
iLanger & Willacv I l(200ll) have presented a n estimate of 
the co re density profile from FIR data and iLevin et al. I 
l|200l[) have estimated an upper limit to its magnetic field 
of 100 IJ,G from Zeeman observations. Tafalla et al. (2002, 
TMCWC hereafter) presented line and continuum obser- 
vations of L1498, L1517B, and three other dense cores, 
and found a systematic pattern of chemical differentia- 
tion: species like CS and CO are almost absent from the 
core interiors while NH3 and N2H+ are present even in 
the densest parts. In this paper, we improve the study of 
L1498 and L1517B using new multiline, high angular res- 
olution observations and creating for each core a model 
of the internal conditions that simultaneously explains all 
our Ci^O, C^^O, CS, C^^S, N2H+, and NH3 data. As a 
result, we confirm the presence of a systematic abundance 
pattern, and show that the two cores are very close to 
isothermal spheres both in their density profiles and physi- 
cal parameters (temperature and turbulence). At the same 
time, we infer that the core contraction history has been 
asymmetrical, both from the internal kinematics and the 
chemical composition. In a follow-up paper (Tafalla et al. 
2004, in preparation) we will present a full chemical survey 
of these two cores. 



2. Observations 

We observed L1498 and L1517B in the 1.2 mm contin- 
uum with the IRAM 30m telescope in 1999 December and 
2000 December. We used the MAMBO 1 bolometer array 
ijKrevsa et al. lll'998(l in on-the-fly mode with a scanning 
speed of 4" s~^, a wobbler period of 0.5 s, and a wobbler 
throw of 53". The data were corrected for atmospheric 
extinction using sky dips done before and after each in- 
divid ual map, and they were reduced using the NIC soft- 
ware (jBroeuiere et al. lll996l) . A global cahbration factor 
of 15000 counts per Jy beam~^ was estimated for both 
1999 and 2000 epochs from observations of Uranus. The 
bolometer central frequency is 240 GHz and the telescope 
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beam size is approxima tely 11". (The 1999 data have al- 
ready been presented in lTafalla et al. l l2002.) 

We observed L1498 and L1517B in Ci*^O(l-0), 
Ci«0(2-1), Ci^0(2-1), CS(2-1), CS(3~2), &^S{2-1), 
C34s(3-2), and N2H+(l-0) with the IRAM 30m tele- 
scope in 2000 October. L1498 was mapped in all Hues and 
L1517B was mapped in aU but 0^48(2-1) and 0^48(3-2), 
which were so weak that only the central position was ob- 
served. Additional N2H+(3-2) observations were carried 
out in 2001 April. All data were taken in frequency switch- 
ing mode with frequency throws of approximately 7 or 14 
MHz. As a backend we used the facility autocorrelator, 
split into several windows to achieve velocity resolutions 
between 0.03 and 0.06 km s~^. The telescope pointing was 
corrected using frequent cross scans of continuum sources, 
and the total errors were found to be smaller than 4" 
(rms) . The data were calibrated using the standard chop- 
per wheel method and converted to the main beam bright- 
ness temperature scale. The full width at half maximum 
(FWHM) of the telescope beam is inversely proportional 
to the frequency, and ranges from approximately 26" at 
the (lowest) frequency of N2H+(l-0) to 11" at the (high- 
est) frequency of N2H+(3-2). 

As stressed in previous studies I Lee et ai~1 l200lt 
iTafalla et al. I l20o3) . the modeling and interpretation of 
the very narrow spectral lines from starless cores like 
L1498 and L1517B (FWHM « 0.2 km s^^ w 70 kHz at 
100 GHz) requires using rest frequency values more ac- 
curate than those in standard compilations. Fortunately, 
current laboratory work is starting to produce frequency 
determinations with uncertainties of the order of 1 kHz, 
and whenever such a measurement exists we have pre- 
ferred it over any other value, including our ow n previous 
astronomical determinations l|Lee et al. II2OOII) . The rest 
frequencies used in this paper are indicated in Table 1. 
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Fig. 1. 1.2 mm continuum maps of L1498 (left) and L1517B (right). Note the central concentration and symmetry of 
the emission. In order to enhance the sensitivity, each map has been convolved with a 20" gaussian. First contour and 
contour spacing are 2 mJy/ll"-beam, and central coordinates are (aigso — 4'^7™50!0, Jigso — 25°02'13'.'0) for L1498 
and (ai95o = 4'^52'"7!2, Jigjo = 30°33'18'.'0) for L1517B. 



3. Results 

3.1. Continuum data and density profiles 

In contrast with the molecules, whose emissivity is highly 
sensitive to excitation and chemistry (see below), the dust 
in a starless core is expected to have an emissivity approx- 
imately constant with radius (see also below for caveats). 
This property allows a rather straightforward conversion 
of the o bserved continuum emission into a density pro- 
file fe.g..lWard-Thompson et al. llOoSlAndre et al. Il996i- 
lEvans etaLTl200l[) . and for this reason, we start our anal- 
ysis of L1498 and L1517B by studying their dust emission. 

Figure 1 presents 1.2 mm continuum maps of L1498 
and L1517B, and shows that in both cores the mm emis- 
sion is centrally concentrated. In L1498, the emission ap- 
pears slightly elongated SE-NW, while L1517B is rounder 
with a weak extension to the SW (and a weaker one to 
the NE). Neither of the cores presents evidence for an un- 
resolved point-like component, indicative of an embedded 
object, so they appear as bona fide starless cores. 

To estimate the density distributions of L1498 and 
L1517B from the data shown in Fig. 1, we follow the same 
procedure as in TMCWC (where the December 1999 part 
of this dataset was presented). Thus, we approximate each 
core as a spherical system, and derive a continuum ra- 
dial profile by taking an azimuthal average of its map. As 
L1498 is slightly elongated, we average the data of this 
source in ellipses of position angle —40° and aspect ratio 
b/a = 0.6, while for the rounder L1517B, we average the 
data in circles. These radial profiles are presented in Fig. 
2 both in linear and log-log scales. (Core centers are as- 
sumed to be at (10", 0) and (-15", -15") for L1498 and 
L1517B, respectively.) 



To model each continuum radial profile we need to as- 
sume a dust temperature law. Theoretical models of the 
dust temperature in cores predict a slight i nward decrease 
(iMathis et al. I Il983l: lEvans et al.1 lioOlt IZucconi et al. I 
l200l|) . whichTor the case of the L1498 and L1517B cores 
is estimated to range between 11 K in the outer layers 
and 8 K in the innermost part ( Evans et al. lliooH their 
model for L1512). Here, for simplicity (and from our NH3 
analysis in Section 3.3), we will assume that both cores 
have a constant dust temperature of 10 K, and we will 
use this value to estimate their density radial profiles (see 
alsoE^anecr & Willacv 200^. 

In addition, we assume a constant dust emissivity of 
Ki.2mm = 0.005 cm^ g~ \ for consistency with previous 
work (e.g. jAndre et al. l liggG). Following TMCWC, we 
fit the data using a density profile of the form 



n(r) = 



no 



l + (r/ro)"' 



where no is the central density, rg is the radius of the inner 
"flat" region (2ro is the full width at half maximum), and 
a is the asymptotic power index. For each choice of no, ro, 
and a, we calculate the emergent continuum distribution 
by integrating the equation of radiative transfer in the 
optically thin case, and simulate an on-the-fiy bolometer 
observation of such a distribution using the same instru- 
mental parameters as those used at the telescope. The 
radial profile of this simulated observation is compared 
with the measured radial profile, and the free parameters 
no, ro, and a are modified until model and data agree. As 
shown by the solid lines in Fig. 2, the final choice of free 
parameters (Table 2) gives very reasonable fits for both 
cores, and indicates that L1517B is twice as dense at its 
center than L1498, while it has a FWHM approximately 
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Table 2. Density fits from continuum data 
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Fig. 2. Radial profiles of 1.2 mm continuum emission 
(from the maps in Fig. 1). The squares represent observed 
data and the lines are the predictions from two density 
models: the solid line is a simple analytic model and the 
dashed line is an isothermal model (see text) . Both linear- 
linear and log-log plots are shown. 

half that of L1498 (so both cores have similar central col- 
umn densities: 3-4x10^^ cm~^). We note that the fit for 
L1517B is the same as derived in TMCWC, while the fit 
for L1498 is only slightly different (by less than 10% for 
radii smaller than 110" and by less than 30% in the range 
110-150"). 

Although our analytical density distributions accu- 
rately fit the continuum data, their motivation is purely 
phenomcnological. As detailed in Appendix A, however, a 
particular choice of parameters in this density family pro- 
vides an excellent approximation to an isothermal func- 
tion (or Bonnor-Ebert sphere), and this choice happens 
to be the best fit for L1517B (see Appendix A for fur- 
ther details). This "coincidence" motivated us to explore 
isothermal fits to the L1517B and L1498 radial profiles 
as an alternative to the analytical models presented be- 
fore. For this, we have used the numerical solution to th e 
isothermal function by IChandrasekhar fc Wares! l|l949l|) , 
and searched for the best fit to each core by varying the 
central density and the isothermal sound speed as inde- 
pendent free parameters. In this procedure, we have fol- 
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lowed the steps explained in the fitting with analytical 
models including the on-thc-fly simulation, and the best 
fit parameters are given in Table 2. As can be seen in Fig. 
2, the isothermal models (dashed lines) provide reason- 
able fits to the observed continuum data, especially for 
the rounder L1517B, where the isothermal and analytic 
models are indistinguishable (see also Appendix A). For 
the more elongated L1498 core, the isothermal model does 
not fit the radial profile as well as the analytic function, 
but its deviation from the data is never very large (~ 25% 
at the center), and is less than the deviation of this core 
from spherical geometry (about 40%). Our fits of L1498 
and L1517B, therefore, add to recent evidence that the 
density profiles of starless cores can be modeled exactly 
or approximately by isothermal spheres jBacmann et al. I 
l200ntlAlves et al. "200 1^ 'Evans et al. "20Ql|) ~ 

As shown by Ebert (1955) and Bon nor I l|l956() . a self 
gravitating isothermal sphere is stable only if it is bounded 
from the outside and its center-to-edge density contrast 
does not exceed a limiting value close to 14, or equiv- 
alently, its dimensionless radius (5) does not exceed the 
critical value 6.45. Although neither L1498 nor L1517B 
show a clear boundary at which the core merges with an 
external medium, we can assess possible core stability by 
calculating the critical radius for each core. For the less 
centrally concentrated L1498 core, the critical radius oc- 
curs at 150", which is at the very edge of our map (defined 
by the sensitivity of our continuum measurements). For 
L1517B, the critical radius occurs at 100"; so if the isother- 
mal core continues to the edge of our map, the core would 
seem unstable. This situation is similar to that found by 
Bacmann et al. ( 2000) in their study of a sample of star- 
less cores with ISOCAM. We defer further discussion of 
the isothermal model and its stability to section 4.1. 

3.2. Overview of the molecular line data and model 
parameters 

To derive the temperature, kinematics, and chemical 
structure of L1498 and L1517B, we now analyze their 
molecular emission. Before starting the detailed radiative 
transfer modeling, we first present a comparison between 
the integrated maps of all species and identify the most 
general trends. As Figure 3 shows (see also TMCWC), the 
maps of each core display a striking dichotomy — almost 
to the point of anti correlation — between the distribution 
of the different tracers: the 1.2 mm continuum, N2II+, and 
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Fig. 3. Integrated intensity maps for continuum and lines observed toward L1498 and L1517B. Note the systematic 
dichotomy between centrally-peaked and centrally-depressed morphologies. For each core, the top row shows the 
maps of centrally-peaked tracers, 1.2 mm continuum, N2H+, and NH3, and the two lower rows present the maps of 
centrally-depressed species, isotopomcrs of CO (middle row) and isotopomers of CS (bottom row). In each map, the 
lowest contour and the contour interval are equal, and for each line, the same contour choice has been used for both 
L1498 and L1517B. Lowest contours are (in K km s'^): 0.15 for 0^^0(1-0) and (2-1), 0.05 for 0^^0(2-1), 0.25 for 
N2H+(1 0), 1.5 for NH3(1,1), 0.15 for CS(2-1), 0.1 for CS(3-2), and 0.05 for C^'^S{2 1). Continuum maps and central 
positions as in Fig. 1. The C^^0(2-l) data have been convolved with a 35" (FWHM) gaussian to improve S/N. 
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NH3 emission on the one hand are compact and centrally 
concentrated, while the CO and CS isotopomer emission 
on the other hand are extended and almost ring-like. In 
the larger L1498 core, the CO/CS ring is well resolved 
spatially, and it appears fragmented with a bright peak 
toward the SE. In the more compact L1517B core, the CO 
isotopomers form a well-defined half ring toward the west, 
while the CS emission is more concentrated, although it is 
still offset toward the west from the 1.2 mm/N2H+ maxi- 
mum. 

As modeled in detail in the next sections, the ring dis- 
tributions of the CO and CS isotopomers result from a 
sharp drop in the abundances of these molecules toward 
the dense centers of the cores, most likely caused by the 
freeze out of th ese molecules onto dus t grains (for addi- 
tional cases, seelkuiper et al. I»1996; .Willacy et al. lllQQSt 



tionai cases, seelKuiper et ai. Itiuijo: ■Willacy 
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lBacma,nn et al. Il2nn2l: TMCWC). N2H+ and NH3 
on the other hand, seem immune to this process at typical 
core densities and they are present in the gas phase all 
throughout the core. 

The goal of the following sections is to derive self 
consistent models of the internal structure of L1498 and 
L1517B that fit simultaneously all line observations. For 
this, we model the cores as spherically symmetric sys- 
tems and predict the emission from each molecule with 
a Monte Carlo radiative transfer code based on that of 
iBernes I l|l979l) . The physical properties of these models 
are fully defined by radial profiles of density, tempera- 
ture, turbulent velocity, and systematic velocity, and the 
same parameters are used to model all molecular lines. 
As core density profiles, we use the analytic functions de- 
rived from the continuum data, and the rest of the gas 
parameters are derived from the molecular line data. In 
our derivation of these parameters we have favored sim- 
plicity over detailed fine tuning, and here we summarize 
the main results (see following sections for more details 
and Table 3 for a parameter list). To derive gas temper- 
ature profiles, we have used the NH3 data, and we have 
found that this parameter is constant with radius in both 
cores (section 3.3). The turbulent profiles are also derived 
from the NH3 data, and they are constant with radius, 
as indicated by the radial profiles of linewidth (section 
3.7.1). The profiles of systematic velocity have been de- 
rived from the combined fit to all line shapes, as different 
species trace the gas velocity at different radii depending 
on optical depth. For both L1498 and L1517B, we find 
that the velocity field in the core central part is constant, 
as indicated by the NII3 data, and that there is a slight 
gradient toward the outside traced by the velocity of the 
CS self absorptions (section 3.6). For simplicity we have 
modeled this outer gradient with a linear function. 

When modeling the self absorbed lines of CS (section 
3.6), HCO+, and HON (Tafalla et al. 2004 in preparation), 
the shape of the profiles depends very sensitively on the 
gas conditions at very large radn (> 4 x lO" cm « 190"). 
These conditions, however, cannot be safely extrapolated 



from our dust continuum or NII3 data, as they correspond 
more to the ambient cloud or envelope than to the dense 
core itself, and favoring again simplicity, we have mod- 
eled them using constant gas parameters determined from 
an approximate fit to the line self absorptions of different 
molecules. For both L1498 and L1517B, we have intro- 
duced a discontinuity in the density profile at 4 x 10^^ cm, 
where both cores reach a density of approximately 3 x 10'^ 
cm~^, and we have extended these profiles with envelopes 
of constant density equal to lO'^ cm~'^ up to a radius of 
8 X 10^^ cm. To match the self absorption features, we 
have set the velocity of this envelope equal to zero for 
L1517B and to 0.35 km for L1498 (inward motion), 
and the turbulent component equal to 4 and 5 times the 
central turbulence, for L1517B and L1498 respectively. As 
we will see below, the inward-moving envelope component 
in L1498 seems only present in the front part of the cloud 
and most likely represents an unrelated foreground cloud, 
so we have neglected the contribution of the back side of 
the envelope in this object. It should be stressed that our 
parameterization of the cloud envelopes is highly simpli- 
fied and it has been carried out with the only purpose of 
fitting the shape of the self absorbed lines. These envelopes 
have no effect on the emission of the thin lines used to de- 
termine the abundance profiles of each species, so their 
role in the modeling is almost cosmetic. Detailed mapping 
of the environment surrounding the cores is necessary to 
further constrain the envelope properties. 

Having fixed the physical parameters of the cores 
(summarized in Table 3), the only free parameter left to 
fit the emission of each species is the radial profile of its 
abundance. To determine this profile, we have searched 
for Monte Carlo solutions of the radiative transfer that 
fit simultaneously the radial profile of integrated inten- 
sity and the central spectrum of each observed transition. 
This approach is similar to that of TMCWC, with the nov- 
elty that the present analysis uses higher spatial resolution 
observations, and that for each species we simultaneously 
model multiple transitions. This new analysis results in a 
more detailed and accurate derivation of the physical and 
chemical structure of the cores. 



3.3. NH^: constant temperature and central 
abundance enhancement 

TMCWC used NH3 data to determine the gas kinetic tem- 
perature of several cores including L1498 and L1517B. 
Given the importance of this parameter and the re- 
cent calculation s of the expected temperature structure 
of dense c ores fe vans et al. 1 120011: IZucconi et al. I l200lt 
ICalli et"al. . 2002) , here we present an improved NH3 anal- 
ysis using the new density profiles of section 3.1 and a 
more accurate temperature determination method. Our 
new met hod is an updated version of th e classical anal- 
ysis (e.g.. IWalmslev fc Ungerechts Ill98:^ calibrated with 
Monte Carlo radiative transfer models, and as detailed in 
Appendix B, it provides accurate temperature determina- 
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Table 3. Monte Carlo model fits 



Parameter 


L1498 


L1517B 


Note 


Core 


Tmax (cm) 


4 X 10^^ 


4 X 10" 




Vlsr (km s"^) 


7.8 


5.8 




dV/dr\in (km s~^ pc~^) 








(1) 


dV/dr\out (km s~^ pc"^) 


-0.5 


2.0 


(1) 


Tbr (cm) 


f.75 X fO^^ 


1.5 X 10" 


(1) 


AVnt (kms-i) 


0.125 


0.125 




n (K) 


10 


9.5 


(2) 


Envelope 


Tniax (cm) 


8 X 10" 


8 X 10" 




Vlsr (km s~^) 


7.45 
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tions for realistic core conditions. With its help, we first 
estimate the L1498 and L1517B temperature profiles di- 
rectly from the NH3 data, and then use the full Monte 
Carlo model to fit the observed radial profiles and spectra 
deriving NH3 abundance profiles (and demonstrating the 
self consistency of the temperature determination) . 

To derive the temperature profiles of L1498 and 
L1517B, we first estimate the T^^ rotation temperature for 
each observed position using the NH3 lines. This involves 
fitting the hyperfine (hf) structure of the NH3(1,1) spec- 
tra (where we find that the hf components are in LTE), 
fitting single gaussians to the NH3(2,2) spectra, and esti- 
mat ing the ratio between the (2,2) and (1,1) populations 
(see lBachiller et al. ll987l for details). Then, we use the an- 
alytical expression derived in Appendix B to convert each 
rotation temperature into a gas kinetic temperature. The 
resulting temperature radial profiles, presented in Fig. 4, 
show that in each core the gas kinetic temperature is to 
very good approximation constant with radius and almost 
equal to 10 K with an rms of about 1 K. 

The temperature profiles in Fig. 4 represent line-of- 
sight averages spatially smoothed with a 40" gaussian 
(telescope resolution). To test the possibility that a real- 
istic non-constant temperature distribution could mimic 
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Fig. 4. Radial profiles of gas kinetic temperature for 
LI 498 and L1517B derived from NH3. The squares rep- 
resent real data and the solid lines represent the result of 
applying the same NH3 analysis used for the data to spec- 
tra generated with a Monte Carlo model having constant 
temperature (10 K for L1498 and 9.5 K for L1517B). The 
dashed lines are the expected NHs-derived Tk profile for 
two temperature equilibrium models (see text for details) . 



the constant temperature profile of our NH3 analysis, we 
have run a Monte Carlo NH3 model fo r L1517B using the 
temperature law recently calculated bv lGalU et all ll2002^ 
(and kindly adapted to L1517B by Daniele Galli). The 
synthetic NH3 spectra from this model have been analyzed 
as the real data and used to compute an NH3-derived tem- 
perature profile. As the upper dashed line in Fig. 4 shows, 
a model with the "s tandard" cosmic r ay ionization rate of 
^ = 3 X 10"" s-i llGoldsmithll200l|) . overestimates the 
core temperature at all radii (in fact, it overestimates the 
(2,2) intensity by more than a factor of 2). However, re- 
ducing C by a factor of 5 (lower dashed curve) produces a 
temperature profile consistent with the observations (such 
a low C, value is also favore d by the chemical m odels of 
the LI 544 core in T aurus bv lCaselh et al. ll2002bl but see 
McCall et al. l200,l for the opposite conclusion in a nearby 
cloud). This low ( temperature profile is in fact almost 
constant, as it has a maximum-to- minimum variation of 
less than 1.5 K in the region of interest. We thus conclude 
that the gas kinetic temperature in the L1498 and L1517B 
cores is constant within about 10%, and that this property 
is consistent a with realistic treatment of the heating and 
cooling of dense gas (assuming a low C value). 

Given the above results, we model both L1498 and 
L1517B as having constant temperature, and for each core, 
we search for the NH3 abundance profile that fits simul- 
taneously the observed (1,1) and (2,2) radial profiles of 
integrated intensity together with the central spectra of 
both transitions. We do this by using our full non LTE 
Monte Carlo radiative transfer analysis, and convolve the 
results with a 40" gaussian to simulate an observation with 
the Effelsberg antenna (see Appendix B for details). For 
L1517B, our best model is identical to that in TMCWC, 
as we use the same density profile. This means that the 
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Fig. 5. Radial profiles of integrated intensity and central spectra of NH3(1,1) and (2,2) towards L1498 and L1517B. 
The squares in the radial profiles and the histograms in the spectra represent real data. The solid lines represent the 
results of Monte Carlo radiative transfer models with constant temperature and a central NH3 abundance enhancement 
(see text). Note how they simultaneously fit the radial profiles and the central spectra (they also fit the NHa-derivcd 
temperature profiles, see previous figure). The dashed lines in the radial profiles are models with constant NH3 
abundance, chosen to fit the observations at large radius (^ 60"). Note how these flatter profiles cannot fit the 
observed central emission. The velocity scale in L1517B has been shifted by 2 km for presentation purposes. (Note 
that the 1-sigma uncertainty in the integrated intensity points is of the order of 0.2 K km s~^ for NH3(1,1) and 0.8 K 
km s~^ for NH3(2,2) multiplied by 20, and it is therefore similar to the size of the squares.) 



gas temperature has a constant value of 9.5 K and that 
the para-NHa abundance increases inwards following the 
law X{r) = Xo{n{r)/nof, where Xq = 1.7 x 10"^ and 
(3=1 (the above expression is only valid for X > 10~^). 
Our best fit model for L1498 is slightly different from that 
in TMCWC because of our improved density profile and 
core peak determination from the new continuum data. 
This model has a constant temperature of 10 K and a 
para-NHs abundance that increases inwards following the 
same law as in L1517B, but this time with Xq = 1.4 x 10~^ 
and /3 = 3, and again with a 10"^ threshold. (Note that 
if the ortho-to-para ratio of NH3 equals 1, the total NH3 
abundance is 2 times the measured para-NH3 abundance.) 
As Fig. 5 shows, these NH3 models fit simultaneously the 
radial profile of integrated intensity and the central emer- 
gent spectra for both NH3(1,1) and (2,2). They also fit 
NH3-derived temperature profiles of Fig. 4 (see solid lines). 

Although initially unexpected, the higher central NH3 
abundance in L1498 and L1517B is a necessary element of 
the radiative transfer models. To illustrate this point, we 
also present in Fig. 5 (dashed lines) the results of constant 
abundance models chosen to fit the emission at interme- 
diate radii (w 70"). As can be seen in the figure, these 
models predict radial profiles that are too flat for both 
NH3(1,1) and (2,2), and are therefore inconsistent with 
the data. The higher central abundance, thus, is unavoid- 
able, and sets apart the NH3 molecule among all other 



observed species (see Tafalla ct al. 2004, in preparation, 
for an extensive chemical survey of L1498 and L1517B). In 
section 4.2 we discuss the possible origin of this behavior. 

3.4. N2H^ : constant abundance 

Like NH3, N2H+ presents centrally peaked distributions 
in our maps of Fig. 3, and this suggests that the N2H+ 
abundance does not vary greatly over the core interiors. 
In this section we quantify this impression by solving the 
N2H"'" radiative transfer with our Monte Carlo code. This 
code now needs to take into account the hyperfine struc- 
ture introduced by the two N atoms, because (in contrast 
with the almost thermalized NH3) the N2H"'" hyperfine 
structure affects the line trapping and the level excita- 
tion. Unfortunately, an exact treatment of the N2H+ ra- 
diation transfer with hf structure is not yet possible for 
two reasons. First, the collision rates between the individ- 
ual hf components are not known, and only approximate 
hf-averaged values are available (see below). Second, the 
hf structure introduces the possibility of partial overlap 
between lines, complicating the problem. 

TMCWC presented a simplified solution of the N2H+ 
transfer which ignored the hf structure in the excitation 
calculation but considered it in the computation of the 
emergent spectrum. As mentioned there, the main prob- 
lem with this approach is that it neglects the decrease 
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Fig. 6. Radial profiles of integrated emission (left) and central emerging spectra (right) for N2H+(l-0) toward L1498 
and L1517B. Observed data are represented by squares in the radial profiles (sum of all hf components) and by 
histograms in the spectra (main beam brightness temperature). The solid lines represent the results from best fit 
Monte Carlo models: for L1517B, the model has a constant N2H+ abundance of 1.5 x 10~^°, and for L1498 the model 
has a constant abundance of 1.7 x 10^^° with a drop of a factor of 3 for radii larger than 1.8 x 10"'^'^ cm. The dashed line 
in the radial profile of L1498 shows the expected emission for a constant abundance model without outer drop. The 
velocity scale in L1517B has been shifted by 2 km s~^ for presentation purposes. (Note that the 1-sigma uncertainty in 
the integrated intensity points is of the order of 0.05 K km s^^ and it is therefore smaller than the size of the squares.) 



in trapping caused by the splitting. Here, to take this ef- 
fect into account, we improve the model and treat explic- 
itly the hf splitting in the low J transitions (up to J=2), 
where the trapping decrease is expected to be important, 
while we assume the rest of transitions are degenerate (up 
to J=6). This improved approach requires increasing the 
number of levels considered in the Monte Carlo code to 23 
and the number of transitions to 44 (no line overlaps are 
consi dered). For col lision rates, we use the HCO+-H2 rates 
from lFlowerl (|l£9?) and derive the individual coefficients 
between N2H"'" hf components by assuming that collisions 
do not distinguish between hf levels. By using these sim- 
plified collision rates and neglecting line overlaps we lose 
the ability to reproduce the non LTE excitation known to 
occur m the J=l-0 line ijCaselh et al. Ill995j) . This exci- 
tation anomaly, however, seems to be a 10% effect in the 
level population, and therefore is of little consequence for 
our abundance estimate. 

To better constrain the N2H+ abundance profiles, we 
fit simultaneously the N2H+(l-0) radial profile of inte- 
grated intensity and the emergent spectra toward the core 
center for both N2H+(l-0) and N2H+(3-2). We start our 
analysis by exploring models with constant abundance, 
which wc find fit the data rather well. As illustrated in 
Fig. 6 and 7, a constant abundance model with X(N2H+) 
= 1.5 x 10"^" fits automatically all L1517B observations 



and no additional change is required in the abundance 
profile to fit this source. (Note that our model spectra are 
slightly brighter than the observations despite having the 
same integrated area because the real N2H+ lines seem to 
have additional broadening, most likely caused by unre- 
solved hyperfine structure, see section 3.7.1). For L1498, a 
constant abundance model with X(N2H+) = 1.7 x 10~^° 
fits well both the J=l-0 and 3-2 spectra (not shown) and 
almost fits the J=l-0 radial profile, although it tends to 
overestimate the intensity at large radius (dashed lines in 
Fig. 6). To improve this fit, we have modified the constant 
abundance model by introducing an abundance drop of a 
factor of 3 for radii larger than 1.8 x 10^''' cm (85" at 140 
pc). This new model fits all L1498 data rather well, as can 
be seen by the solid lines in Fig. 6 and 7 (the slight devi- 
ation from flat top at the center of L1498 probably arises 
from a similar deviation in our density profile, see Fig. 2). 

In contrast with our conclusion that the N2H+ abun- 
dance is constant o r close-to-cons tant at the centers of 
L1498 and L1517B. iBergin et al. I f2002) find a factor of 
2 drop in the N2H+ abundance at the center of the B68 
core. Such a significant drop seems excluded in L1517B 
or L1498, not just from our modeling but from the lack 
of any significant drop in the N2H+ emission near the 
continuum peak (see maps in Fig. 3). This contrast with 
the behavior in B68 is especially significant in the case of 
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Fig. 7. Observed N2H+(3-2) spectra (histograms) and 
predictions (lines) for the same models shown in Figure 
6. The models derived to fit the 1-0 transition also fit the 
3-2 spectra. (Note that the multiple peaks in the spec- 
tra arise from hyperfine structure, and do not represent 
multiple velocity components.) 



L1517B, as both B68 and L1517B have central densities 
of abou t 2 X 10^ cm~'^ an d are well fitted by isothermal 
models llMves et a1. 112001 [ our section 3.1). Given the im- 
portance of N2H"'" as a tracer of dense gas, further work 
should be carried out to understand the discrepancy be- 
tween these cores. 

3.5. CO isotopomers: central depletion 

The central holes in the CO isotopomer maps of Fig. 3 
indicate graphically that the abundance of these species 
drops systematically toward the peak of both L1498 and 
L1517B. To determine quantitatively these abundance 
drops, we use again the Monte Carlo radiative transfer 
model to fit the observations with a variable abundance 
profile (for further details on the Monte Carlo modeling, 
see TMCWC). To better constrain the abundances and to 
provide a self consistency check to our calculations, we use 
simultaneously all the CO isotopomer data we have avail- 
able. Thus, in addition to fitting the 01^0(1-0), C^^0{2- 
1), and 0^^0(2-1) data from the IRAM 30m telescope, 
we use the lower angular resolution C^^O(l-O) data from 
FCRAO already presented in TMCWC. Also, taking into 
account previous work on isotopic abundances in the ISM, 
we assume a co nstant C^*0/C^^O abundance ratio of 3.65 
llPenzias 1 ll98ll) . In this way, our determination of the CO 
isotopomer abundances reduces for each core to finding 
a single abundance profile that fits simultaneously four 



radial profiles of integrated intensity and four emergent 
spectra. 

Before attempting to reproduce the observations, we 
test the conclusion that constant abundance models can- 
not fit the data from neither L1498 nor L1517B. To do 
this, we fix the C^'^O abundance so that the models fit 
the observations at a radius of 100", and compare the 
model predictions with the observations at smaller radii. 
As the dashed lines in the radial profiles of Fig. 8 show, the 
constant abundance models fail to fit the data for all CO 
isotopomers in both L1498 and L1517B. They do so by a 
wide margin (more than a factor of 2), and this proves that 
constant abundance C^^O/C'^^O models are inconsistent 
with observations. A sharp drop of abundance toward the 
center is therefore required. 

To model the CO central abundance drop, we 
have explored different abundance profiles, starting with 
the exponential forms introduced by TMCWC {X ~ 
exp(— 7i(r))). Our higher angular resolution data, however, 
suggest that these forms do not provide a fast enough 
abundance drop to fit the integrated intensity near the 
core centers, and that faster decreases are needed to re- 
produce the data (note that because of the central flatten- 
ing of the density profiles, the exponential forms do not 
necessarily result in dramatic abundance drops). Thus, we 
have explored alternative abundance profiles, like squares 
of the exponential forms or central holes. Unfortunately, if 
the abundance toward the core center is low enough that 
the emission is dominated by the outer layers (order of 
magnitude drop), the exact form of the abundance profile 
is not very critical. Given this degeneracy of the solution, 
we have chosen the simplest form, a (step) central hole, 
and used it as basic model for the CO (and CS) abundance 
profiles in L1498 and L1517B. Although in some cases a 
slightly different abundance profile may produce a better 
fit to the integrated intensity radial profile (as in the case 
of CS in L1517B, see below), we have preferred for con- 
sistency and inter-comparability to fit all the abundance 
drops with central holes. 

Fig 8 shows a comparison between the central hole 
models (solid lines) and the observations for both L1498 
and L1517B. These models are characterized by just two 
parameters (external abundance and hole radius, see Table 
3) and provide a simultaneous fit to the four integrated 
intensity radial profiles and the four central spectra of 
0^*0 and C^'^O. Their abihty to fit the data shows that 
the central abundance drop in each core is sudden and 
deep, and that the data are consistent with both L1498 
and L1517B having no CO molecules in their innermost 
region (r < 70"). Although the radius of the central hole is 
not uniquely determined by the observations (it is slightly 
coupled to the initial abundance estimate), the true values 
cannot be very far from our estimates. Thus, we combine 
the radius of the abundance drop (Table 3) with the den- 
sity profiles (Table 2) in order to estimate that the density 
at which CO molecules disappear from the gas phase is 
7.8 X 10^ cm-3 for L1498 and 2.5 x 10"* cm'^ for L1517B. 
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Fig. 8. Radial profiles of integrated omission (left) and central emerging spectra (right) for C^*O(l-0), C^®0(2-1), 
C^^0(1 0), and C^''0(2-l) toward L1498 and L1517B. Observed data are represented by squares in the radial profiles 
and by histograms in the spectra (main beam brightness temperature). The dashed lines represent constant abundance 
models chosen to fit the emission at radius 100" {X{C^^O)= 0.5 x 10""^ in L1498 and 1.5 x IQ-"^ in L1517B) and, as the 
figure shows, they fail to fit the central emission by a large margin. The solid lines are models with the same abundance 
value for large radii but having a central hole with no molecules (r^o/e = 71" for L1498 and 83" for L1517B). These 
models fit simultaneously all transitions at all radii. A C-'^^O/C-'^'^O isotopic ratio of 3.65 is assumed in both cores. 
(Note that the 1-sigma uncertainty in the integrated intensity points is of the order of 0.03 K km s~^ and it is therefore 
similar to the size of the squares.) 
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3.5.1. Deviations from spherical symmetry 

Our radiative transfer analysis has assumed spherical sym- 
metry, and our Monte Carlo models have been used to fit 
azimuthal averages of the observed emission. This means 
that our radial profiles with a central abundance drop 
represent spherical averages of the distribution of CO 
molecules in each core. As the maps in Fig. 3 show, how- 
ever, the distribution of CO emission in both L1517B and 
L1498 is not fully circular (or elliptical), and superposed 
upon the central emission drop, there is a pattern of sec- 
ondary peaks and valleys. In L1517B, for example, the 
CO emission is brighter toward the west, and in L1498 
there are bright spots toward the southeast and the west. 
The asymmetries in the CO emission appear in all tran- 
sitions but do not have counterparts in the continuum or 
NH3 emission (Fig. 3). This suggests that they do not arise 
from asymmetries in the distribution of density or temper- 
ature, but that they result from true non-azimuthal asym- 
metries in the distribution of CO abundance. The detailed 
modeling of these asymmetries is outside the scope of our 
spherically symmetric radiative transfer analysis, but a 
simple inspection of the maps in Fig. 3 suggests that the 
abundance variations needed to produce them are at the 
level of a factor of two. These variations, therefore, are 
only perturbations on the (order of magnitude) drop in 
abundance at high densities, but as we will see in section 
4.3, they offer interesting clues concerning the process of 
core formation. 

3.6. CS: central depletion 

We finish our abundance analysis studying the distribu- 
tion of the CS molecule. Unfortunately, its emission is 
more difhcult to model because the CS lines from both 
cores are very optically thick and strongly self absorbed. 
This can be seen in Fig. 9, where the main CS isotopomer 
lines appear asymmetric and double-peaked while the 
spectra of the rare C'^^S species are single gaussians. These 
self absorbed spectra depend strongly on the physical pa- 
rameters of the low-excitation outermost layers that give 
rise to the absorption (r > 0.1 pc) and for which we have 
very little information in terms of density and tempera- 
ture. Thus, although the Monte Carlo code can treat accu- 
rately the transfer of radiation under optically thick con- 
ditions, our CS modeling is limited by our ignorance of the 
gas properties in the outer core regions. If these properties 
deviate from the extrapolation of our simple models, the 
radiative transfer solution for the main isotopomer of CS 
is compromised. To minimize this effect, we have comple- 
mented our CS observations with C'^^S data. 

As with CO, we test the CS models simultaneously 
against our full CS isotopomer data set, which consists 
of maps and spectra of CS(2-1) and CS(3-2), and central 
spectra of &^S{2~1) and C^^S{3-2) (plus a fuU C34S(2-1) 
map of L1498). To minimize the number of free parame- 
ters, we set the CS/C^^'S isotopic ratio to its solar value 
of 22.7, since this number seems consistent with ISM de- 



terminations l|Lucas fc Liszt Ill998|) . An acceptable abun- 
dance curve, therefore, must fit simultaneously the 4 emer- 
gent spectra and 2 radial profiles (3 in L1498) shown in 
Figure 9. 

Our first set of Monte Carlo models have spatially con- 
stant CS abundance, set to a value that fits the integrated 
intensity radial profile at a radius of 100". As the dashed 
lines in Fig. 9 show, these constant abundance models pre- 
dict a significant concentration of the emission toward the 
core center, which is not seen in the data. This is more 
evident in L1517B because of its higher central density 
concentration, but is also clear in the flatter and more 
CS-opaque L1498 core, especially in the thinner C^''S(2- 
1) line, where the central deviation between model and 
data is of at least a factor of two. The failure of the con- 
stant abundance models to fit the data proves that the CS 
abundance must drop toward the centers of both L1498 or 
L1517B ( in agreement w ith previous lower resolution re- 
sults, see lKuiner et al. Itl99a and TMCWC). 

To model the central decrease of the CS abundance, 
we have again explored different families of abundance 
profiles. Models with a central hole, for example, seem 
slightly superior to exponential drop models in the case 
of L1498, ahhough both models fit the data in L1517B. 
Unfortunately, a detailed comparison between the differ- 
ent forms of the central abundance drop is difficult because 
of the dependence of the results on the CS self absorption, 
and for this reason, we only present the results of models 
with a central hole. These models provide reasonable fits 
to the observations, as illustrated in Fig. 9, which shows 
that central hole models (solid lines) match simultane- 
ously for each core all radial profiles of integrated intensity 
and all central spectra. 

The simultaneous fit of both the self absorbed CS lines 
and the gaussian C'^^S spectra is a critical test for the mod- 
els. The CS self absorption arises as a natural consequence 
of the combination of high optical depth and subthermal 
excitation in the model low-density outer layers, and this 
makes this feature the most sensitive tracer of the veloc- 
ity field in the outer core (as mentioned in section 3.2). 
The red shifted CS selfabsorption in L1498 indicates a 
contraction motion, while the blue shifted self absorption 
in L1517B req uires expansion (in agreement with p revi- 
ous work, see iKuiper et al. I Il995 iLee et al'~lll999l and 
TMCWC). These velocity gradients, however, do not con- 
tinue all the way to the core center, since otherwise the 
central C'^^S spectra would not have the same linewdith 
(within the noise) as the N2H+ and NH3 lines. Thus, we 
have parametrized the velocity field as constant up to a 
radius of 1.75 x lO" cm (83") in L1498 and 1.5 x lO" 
cm (71") in L1517B, and then as increasing/decreasing 
linearly with radius for larger distances (Table 3). To 
fit the broader CS(2-1) self absorption and the slight 
wing in L1498, we have additionally required that the CS 
abundance increases by a factor of 5 in the outer enve- 
lope (r > 4 X 10^^ cm « 190"). As mentioned in sec- 
tion 3.2, however, this envelope most likely represents a 
foreground red component that also gives rise to a red 
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Fig. 9. Radial profiles of integrated emission (left) and central emerging spectrum (right) for CS(2-1), CS(3-2), 
C^'*S(2-1), and C'^*S(3-2) toward L1498 and L1517B. Observational data are represented by squares in the radial 
profiles and by histograms in the spectra (main beam brightness temperature). The dashed lines represent constant 
abundance models chosen to fit the emission at radius 100" (X(CS)= 3 x 10~^ in both L1498 and L1517B) and, as the 
figure shows, they fail to fit the central emission by a large margin. The solid lines are models with the same abundance 
value for large radii but having a central hole with no molecules (r/jo/e = 48" for L1498 and 55" for L1517B). These 
models fit simultaneously all transitions at all radii. A CS/C'^'^S isotopic ratio of 22.7 is assumed in both cores. (Note 
that the 1-sigma uncertainty in the integrated intensity points is of the order of 0.03 K km s^^ and it is therefore 
similar to the size of the squares.) 



wing in the C^^O spectra (Fig. 8, see also iLemme et all 
Foreground contamination, in fact, may also be the 
cause of the blue shifted self absorption in L1517B, as the 
nearby L1517A core is sligh tly bl ue shifted with respect to 
L1517B (Benson & Mvers1 ll989|) . Thus, the velocity gra- 
dients traced with the CS self absorption most likely re- 
flect the conditions in the region where the cores meet the 
extended ambient cloud (r > 0.1 pc) and not an intrin- 
sic core velocity pattern. Although these motions affect 
strongly the shape of the CS spectra, they seem to involve 
only the core outermost layers and not the dense (probably 
star- forming) part. Their interpretation as global core con- 
traction or expansion patterns seems therefore not war- 
ranted. 



The reasonable fit provided by the central-hole models 
shows that the CS abundance drop in L1498 and L1517B 
is as dramatic as that of CO. This implies that the centers 
of these two cores may lack CS molecules for radii smaller 
than 48" in L1498 and 55" in L1517B (Table 3), and that 
the CS emission in these systems only traces the outer 
core layers. As with CO, we can convert the radius of 
the central hole into a critical density above which CS 
molecules disappear from the gas phase. Using the density 
profiles of Section 3.1, we derive a maximum density of 
8 X 10'' cm-3 for L1498 and 4 x 10'' cm'^ in L1517B, which 
are similar to but slightly larger than those measured for 
the CO isotopomers. Thus, CS seems to freeze-out from 
the gas phase at slightly larger densities than CO, and 
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therefore traces slightly more central parts of cores. The 
innermost core regions, however, are totally invisible in 
CS, even when observed with optically thin C^^'S emission. 

A look at the maps of Fig. 3 shows that the CS emis- 
sion presents the same azimuthal asymmetries found in 
CO: relative maxima toward the SE and W in L1498 and 
toward the W in L1517B. This similarity with the pattern 
of C^^O asymmetries and the fact that in L1498 the asym- 
metry can also be seen in the thinner C'^'*S(2-1) emission 
strongly suggest that these features are real and not the 
result of self absorption effects. Thus, following the argu- 
ments presented in the study of CO (section 3.5.1), we con- 
clude that in each core the distribution of CS abundance 
is not azimuthally symmetric. As in CO, this azimuthal 
asymmetry probably arises from changes in the abundance 
of a factor of a few and is superposed upon the order-of- 
magnitude central abundance drop. A full discussion of 
its possible origin in terms of core formation history is 
deferred to section 4.3. 

3.7. Gas kinematics 

Our modeling of the line profiles toward the center of 
L1498 and L1517B has sampled the gas velocity field along 
the central line of sight of each core. This field has been 
found to have a turbulent component approximately con- 
stant with radius, with a possible increase in the outer 
part of L1498 (but most likely due to another velocity 
component). The line-of-sight systemic velocity seems also 
constant over the center of the cores, although it presents 
gradients in the outer layers: contraction in L1498 and ex- 
pansion in L1517B. In this section, we turn our attention 
to the 2-dimensional (plane of the sky) velocity field as a 
complement to the analysis of the central line-of-sight. To 
do this, we make gaussian fits to the spectra at all posi- 
tions and we study the spatial variations of the fit results. 
As only N2H"*' and NH3 are sensitive to the core dense 
gas, we only use these species in our velocity study. Both 
N2H+ and NH3 have the additional advantage of having 
hyperfine structure, which allows for correction of optical 
depth effects in the linewidth determination. 
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Fig. 10. Radial profiles of intrinsic linewidth for L1498 
and L1517B from NH3 and N2H+ hyperfine-structure fits. 
Note the lack of systematic variations and the extremely 
low dispersion. The long-dashed horizontal lines indicate 
the average value for each core and molecule (very close to 

0. 20 km s""'^ in all cases). The short-dashed lines indicate 
the N2H+ linewidth expected when adding the nonthermal 
component derived from NH3 to a thermal N2H+ compo- 
nent at 10 K (see text for details). To maximize sensitiv- 
ity, only points with integrated intensity larger than 1.3, 

1, 0.6, and 0.6 K km s'^ are presented for L1498 NH3, 
L1517B NH3, L1498 N2H+, and L1517B N2H+, respec- 
tively. (As indicated in the text, the 1-sigma uncertainty 
of the linewidth estimates, both from NH3 and N2H+, is 
of the order of 0.01 km s~^, and it is therefore smaller 
than the size of the squares.) 



3.7.1. Gas turbulence 

We study the spatial distribution of the gas turbulent com- 
ponent using the intrinsic linewidths derived from the hf 
fits and determining their radial profiles as we have done 
before with the integrated line intensities. As the intrin- 
sic linewidth is the sum in qu adrature of a thermal and a 
non thermal component 

(e.g.,|Mmi]Il98l, 

and the ther- 
mal part is constant because the gas kinetic temperature 
is constant (section 3.3), any spatial variation of the in- 
trinsic linewidth has to result from an equivalent spatial 
variation of the non-thermal (turbulent) part. 

Fig. 10 presents the radial profiles of the NH3(1,1) 
and N2H+(l-0) intrinsic hncwidths AV for both L1498 
and L1517B. A simple inspection of these profiles suggests 



that there is no significant trend of increase or decrease of 
linewidth as a function of radius, meaning that any non 
thermal component is also spatially constant over both 
cores. To test for the presence of a hidden slope in the 
linewidth distribution, we have fitted the points in Fig. 
10 with a function of the form AV — a + bR, where R is 
the radius in arcseconds and a and b are two free param- 
eters. In three out of four cases the b coefficient is smaller 
than its rms, hence the fit is consistent with b — and 
the distribution of AV is independent of radius. In the 
fourth case (L1517B N2H+ data), the slope is negative, 
and the linewidth slightly decreases with radius (only by 
30% in 100"). This slight decrease, which can be seen in 
the data at large radii, is most likely due to errors in the 
N2H+ hfs fit (see below) and has no counterpart in the 
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NH3 linewidth distribution. We therefore conclude that 
the non thermal component of the velocity field in both 
L1498 and L1517B is constant with radius for as far as we 
can measure. 

The finding of constant linewidth in L1498 and L1517B 
agrees with the results of iBarranco fc Goodman! lll998l) . 
who found the same behavior in the interior of four 
other dense cores. These authors, however, suggested that 
linewidths may increase at the largest observed radii, a be- 
havior not seen in our data. Even if we include the low sen- 
sitivity points at very large radii, the constant linewidth 
trend of Fig. 10 remains unchanged, and the scatter due 
to noise dominates the plots. We therefore conclude that 
the linewidth does not increase with radius, or if it does, 
it starts increasing at a larger radius than suggested by 
iBarranco fc Goodman I |l998il . 

A puzzling result from our hf analysis is that the 
N2H+(l-0) and NH3(1,1) linewidths are almost equal 
(very close to 0.20 km s^^). This is unexpected because 
the N2H+ molecule is 1.7 times heavier than NH3 and 
thus its thermal component is 1.3 times smaller. If both 
molecules are tracing the same gas, their turbulent compo- 
nents should be equal, and the lighter NH3 should always 
present broader lines. To quantify this problem, we have 
assumed that the NH3 hf fit returns the true linewidth 
(thermal plus turbulent), and we have derived its non- 
thermal part by subtracting a thermal component at 10 K. 
When this nonthermal component is added to the N2H+ 
thermal part, we find that the N2H+ linewidths derived 
from the hf fit are approximately 20% larger than ex- 
pected, as illustrated in Fig. 10, where the dotted line 
is the expected N2H"'" linewidth. 

The origin of the larger N2H+ non-thermal compo- 
nent is unclear. Different linewidths for ions and neutrals 
are expected in the presence of magnetic fields, but the 
observed trend in both L1498 and L1517B (larger non- 
thermal component in the ion line) is opp osite to what is 
expected theoretically ijHoude et al. 112000(1 . A more likely 
explanation is that the N2H"'" hf fit overestimates the 
linewidth by about 20%. This is so because the hf analy- 
sis of both molecules assumes LTE among hf components, 
but this is only correct for the NH3 lines. N2H+ is notori- 
ous for showing non- LTE excitation among its components 
llCaselh et al. 1119951) . and this violates a main assumption 
of the anal ysis. In addition, as o ur data show and was first 
noticed bv lCaselli et al. there is a d screpancy be- 

tween the N2II"'" linewidth derived from the hf fit and the 
linewidth measured by a Gaussian fit to the thinnest (red- 
dest) component. This component appears slightly (5-10 
%) narrower than the hf-fit linewidth, contrary to what 
would be expected if the hf fit truly corrected for optical 
depth broadening. This suggests that N2H+ linewidths de- 
rived from standard hf fits may be slightly overestimated 
(20%), either because of a failure in the LTE analysis or 
because of the presence of additional hypcrfinc splitting 
due to magnetic spin-rotation interaction of the H atom 
(the latter option seems less likely according to recent esti- 
mates of the splitting, Luca Dore private communication). 



Of course, this suggestion does not contradict our conclu- 
sion that the linewidth is constant in L1498 and L1517B, 
as it is confirmed independently (with even less scatter) by 
the more reliable NII3 analysis (Fig. 10). According to this 
analysis, the mean intrinsic (NH3) linewidth is 0.204 km 
s"^ in L1498 and 0.196 km s'^ in L1517. Subtracting ther- 
mal components of 10 K for LI 498 and 9.5 K for L1517B 
(section 3.3), we estimate non thermal turbulent (FWHM) 
linewidths of 0.121 km s'^ in L1498 and 0.113 km s^^ in 
L1517B (note that our Monte Carlo models use 0.125 km 
s^^ as a compromise to fit both NII3 and the broader 
N2H+). 

Even if the turbulent linewidth is constant with ra- 
dius on average, it may still present random variations 
across each core, and to explore this possibility we com- 
pare the measured rms{AV) of N2H+ and NH3 (w lO"'^ 
km s~^) with the value expected from the noise level in 
the data. To calculate this expected rms, we have cre- 
ated at each L1498 and L1517B position a pair of simu- 
lated N2II+ and NH3 spectra with the correct intensity 
and noise level, but with all positions having the same 
linewidth. These simulated data result from averaging all 
spectra of each molecule and core, creating a high S/N 
spectrum which has then been scaled for each position to 
match the observed integrated intensity, with noise added 
to match the observed noise level (the averaging process 
broadens the line due to velocity fields, but by less than 
10%). This artificial data set has then been subject to the 
same hf analysis as the real data, so its rms{AV) provides 
a direct measure of the effect of noise on the measured 
linewidth. Comparing this rms{AV) with that measured 
from the real data, we find that rms{AV) of the real data 
is 1.1 ± 0.3, times the rms{AV) of the simulated data, 
which means that the observed dispersion is consistent 
with noise, and that there is no evidence for an additional 
source of linewidth variations. This result, which agrees 
with the line-of-sight analysis of the Monte Carlo model- 
ing, indicates that the data are consistent with a single 
AV value in the inner core (r < 0.05 pc) within better 
than 0.01 km s^^ (5%). 



3.7.2. Velocity gradients 

The other kinematic parameter derived from the hf fit of 
the NII3 and N2II+ spectra is the line center velocity. As 
with the linewidth, we have made radial profiles of the 
line center velocity for both NII3 and N2H+ in L1498 and 
L1517B, and they are presented in Fig. 11. These radial 
profiles are flat on average, and their scatter around the 
mean is larger than the scatter in linewidth by a about a 
factor of 2. This larger scatter of the line center velocity 
is significant according to the model of simulated spec- 
tra used in the previous section, because the hf analysis 
is more sensitive to the line center velocity than to the 
linewidth, especially for N2H+. Comparing the data with 
the model, we estimate that the observed scatter in the 
line center velocity is 0.03 km s~^, or 5.4±2.5 times larger 
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Fig. 11. Radial profiles of line center velocity for L1498 
and L1517B derived from NH3 and N2H+ hyperfine- 
structure fits. The dispersion of the points indicates the 
presence of bulk internal motions. Intensity threshold as in 
Fig. 10. (As indicated in the text, the 1-sigma uncertainty 
of the velocity estimates, both from NH3 and N2H+, is of 
the order of 0.005 km s~^, and it is therefore smaller than 
the size of the squares.) 
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Fig. 12. Gradients of Vlsr from N2H+ hf fits. The ar- 
rows indicate both the direction (from low to high Vlsr) 
and intensity of the gradient (the arrows in the bottom 
left corners represent gradients of 3 km pc~^). The 
contours are N2H^ integrated intensity maps. 



than it should be if all positions had exactly the same ve- 
locity. This implies that the Vlsr. variations seen in the 
radial profiles of Fig. 11 correspond to real changes in the 
gas velocity across the core and not just to noise. 

To study the origin of the velocity changes, we first 
determine their spatial distribution. Using a program sim- 



ilar to that of ICaselli et al. ] ll2002al) . we have calculated 
at each map position the local velocity gradient, both in 
modulus and direction, by comparing the Vlsr velocity of 
each point with that of its neighbors. A graphical represen- 
tation of the N2H"'" results is shown in Fig. 12, and similar 
plots, but with fewer points, were derived for NH3 (their 
similarity confirms that the velocity fluctuations repre- 
sent real velocity variations). As these maps show, both 
the direction and modulus of the velocity gradient change 
across both L1498 and L1517B, indicating that the gradi- 
ents do not originate from simple global velocity patterns 
like rotation. This lack of global patterns, however, does 
not imply that the velocity variations are completely ran- 
dom, because the gradients are correlated over areas larger 
than the telescope beam (26") or the area used to measure 
individual velocity gradients (40"). They therefore have to 
arise from bulk gas motions that affect sizeable fractions 
of the core material. 

As a first estimate of the importance of these motions, 
we calculate the average velocity gradient for each core. 
Both N2H+ and NH3 give consistent results, especially 
for the position angle (-20 and -75 degrees for L1498 and 
L1517B, respectively). For the modulus of the gradient, 
the NH3 data give a shghtly lower value (10% for L1498 
and 40% for L1517B), so we average the N2H+ and NH3 
values and derive 0.75 and 1.1 km s^^ pc^^ for L1498 and 
L1517B, respectively. Thes e valu es are similar to those 
found by iGoodman et al. I l)l993() in other cores, again 
an indication that L1498 and L1517B are representative 
Taurus cores. Our value for L1498, in addi tion, is in ex- 
cellent agreement with a similar estimate bv iCaselli et al. I 
( {20112 j) . The gradients measured in this manner, how- 
ever, are strict lower limits, as in addition to the average 
across the sky required to estimate the mean gradient, 
the velocity estimate at each position already represents 
an average along the line of sight. This average neces- 
sarily cancels spatial variations of the velocity, given the 
lack of coherence over scales of a core diameter (Fig. 12). 
Unfortunately, this effect cannot be corrected without a 
full model of the core kinematics, and as a simple alterna- 
tive for estimating the velocity differences between parts 
of the core, we measure the spread of Vlsr in Fig- Hj 
which is approximately 0.1 km {= 3x rms). The noise 
contribution to this spread is negligible according to our 
synthetic spectra model (3 x rms — 0.013 km s~^), so 
we can consider that the scatter represents true velocity 
variations. It is possible that part of this dispersion arises 
from a global component (e.g., rotation), but the results 
in Fig. 12 show that this component cannot be dominant. 
Crudely assuming equipartition, we estimate the internal 
core motions to be at least 0.05 km s""'^, although they 
can be as high as 0.1 km s"-'^. 

Although the gas motions are small and subsonic 
(sound speed is 0.19 km s^^ for 10 K), and therefore have 
little effect on the core energy budget, their presence in- 
dicates that the cores have not yet attained a state of 
perfect hydrostatic equilibrium. How these motions origi- 
nate is unclear, but several indications suggest that they 
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Fig. 13. Comparison of "high velocity" N2H"'" emission 
to C'^'^S or CS emission showing their similar distribution. 
This is an indication that the gas velocity structure is re- 
lated to the deviations from spherical symmetry in the 
pattern of molecular freeze out, probably due to asymme- 
tries in the core contraction history (see text). The solid 
lines in the left panels are N2H+ integrated intensity maps 
for velocity intervals 0.2 km s~^ (one linewidth, see Fig. 
10) away from the line center velocity (towards blue for 
L1517B and towards red for L1498) and 0.2 km s~^ wide. 
The C'^'^S or CS maps in the right panels and the full 
integrated N2H+ map (dashed) are as in Fig. 3. 

are related to the process of core formation. First, the 
motions are similar in magnitude to the inwa rd motions 
observed in other cores IILee et al. 1 ll99<il200lh . and these 
are believed to be associated with core contraction (note 
that L1498 has also infall asymmetry, while L1517B has 
the opposite). Second, the time scale of these motions is 
of the order of 1 Myr (0.05 k m s~^ over 0.05 pc o r 75"), a 
typical starless core lifetime (jLee fc Mvers Il999|) . Finally, 
and more important, there is a correlation between some 
velocity components and the non azimuthal asymmetries 
in the abundance of the depleted molecules (CS and CO). 
These abundance asymmetries are most likely related to 
different contraction times of different parts of the core 
(section 4.3), and thus their correlation with the velocity 
pattern suggests that the pattern is related to contrac- 
tion. To show this, we present in Figure 13 (left panels) 
the distribution of "high velocity" N2H+ emission (blue 
for L1498 and red for L1517B) superposed upon the in- 
tegrated N2H"'" emission (see caption for details). As the 
right panels in the figure show, these "high velocity" N2H+ 
components have distributions remarkably close to the dis- 
tribution of bright CS/C'^^S spots (more so than to CO 
probably because CS and N2H+ are sensitive to similar 



densities), which we have seen most likely arise from differ- 
ences in the amount of CS freeze-out. The similar distribu- 
tions are more remarkable because CS and CO hardly co- 
exist with N2H"'", as they deplete at radii for which N2H+ 
is non-detectable, and suggest a correlation of high CS 
abundance with the "high velocity" components. As we 
will see below, these high CS-abundance spots most likely 
arise from material that has been recently acquired by the 
core and has therefore been at high density for shorter 
time th an the rest (i.e., it is "younger" and less depleted, 
see also lKuiper et al.1ll996ll . A peculiar velocity in this 
material is exactly what would have been expected if this 
interpretation is correct. 

4. Discussion 

4.1. Core structure and equilibrium state 

We have seen that the systematic internal motions in both 
L1498 and L1517B are of the order of 0.1 km s~^ and 
therefore subsonic for gas at 10 K (sound speed is 0.19 
km s~^). This suggests that both cores are close to a state 
of hydrostatic equilibrium in which the core self gravity is 
balanced by internal forces. Given the core parameters we 
have derived in previous sections, we now analyze whether 
such an equilibrium state is possible and whether it is 
consistent with the observed core structure. 

The simplest support component in each core is the 
thermal pressure Pt = nkT, where n is the gas density, 
k is the Boltzmann constant, and T the gas kinetic tem- 
perature. This pressure may be supplemented by a non- 
thermal component Pnt — rana'fjrp if the non-thermal 
linewidth arises from turbulent supporting motions (m 
is the mean particle mass and (Tnt is the non thermal 
velocity dispersion, related to the non thermal FWHM 
by AViVT = V8ln2 unt)- The relative importance of 
these two supporting components is given by the ratio 
Pnt / Pt = o'nt / {^"^ / ''^) J which can be easily estimated 
from the core parameters derived in previous sections. 
As both aNT and T are constant over the central 0.1 
pc of L1498 and L1517B, the pressure ratio is constant, 
and has a value of 0.07 for AVpj-T = 0.12 km s^^ and 
T = 10 K (Sections 3.3 and 3.7.1). This low ratio indi- 
cates that the non thermal motions contribute negligibly 
to t he core support over a t least the central 0.1 pc (see 
also iFuller fc Mverslll993l) . and that in the absence of 
support by an ordered magnetic field, thermal pressure is 
the main supporting agent in both L1498 and L1517B (a 
similar conc lusion for the cas e of B68 has been recently 
obtained bv lLada el aTlbOOSl) . 

If thermal pressure dominates core support and the gas 
temperature is constant, the equation of core hydrostatic 
equilibrium reduces to the isothermal c ase, which has a 
classical solution for spherical geometry ijChandrasekhar I 
.1939.) . The addition of a constant non-thermal contribu- 
tion implies only a simple generalization of the isother- 
mal case and does no t change the solution siRiiificantly 
(e.g.. Appendix A of iMcLaughlin fc Pudritz . .1996.) . As 
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we have seen in section 3.1, the density profiles derived 
from the dust continuum emission do in fact have shapes 
very close to those expected for isothermal spheres, es- 
pecially in the case of L1517B. To check whether these 
profiles are consistent with equilibrium between self grav- 
ity and internal (thermal plus non thermal) pressure, we 
can compare the measured velocity dispersion a (=0.185 
km s^^ for gas at 9.5 K and AVnt = 0.11 km s^^) with 
the velocity dispersion required by the isothermal fit to 
the dust emission (section 3.1). Unfortunately, the fit only 
constrains the combined parameter ay^ k^, B{Td), making 
impossible a full check without independent estimates of 
the dust emissivity k^, and temperature Td- ^ To take into 
account this dependence, we re-write the velocity disper- 
sions required from the continuum fits in Table 1 (where 
we assumed = 0.005 cm^ g~^ and Td = 10 K) as 



0.27 (k,/0.005)"0-5 {B{Td)/B{W)) 



-0.5 



km 



for 



L1517B and 0.32 (k^/0.005)-° '5 {B{Td)/B{10))-^-^ km 
for L1498. These values show that if = 0.005 cm^ 
g^^ and Td = 10 K, there is either an extra pressure con- 
tribution which does not affect the observed linewidth or 
the cores should be contracting under their own gravity. 

Although a state of slow contraction is likely in both 
L1498 and L1517B (section 3.7.2), we cannot rule out that 
at least part of the velocity dispersion mismatch refiects a 
wrong choice of dust parameters. For example, the value 
of Td — lOK, chosen so it equals the well-measured gas 
temperature, may in fact be a slight overestimate. This 
is suggested by the models of Galli ct al. (20q3) (their 
Fig. 3), which show that although dust and gas are ther- 
mally coupled at core densities the dust is cooler by about 
2 K. A lower dust temperature requires a larger a, and this 
makes the pressure unbalance more severe, although only 
by a small margin: a dust temperature of 8 K requires a 
1.2 increase in a. More critical in the balance argument 
is the choice of the dust emissivity this parameter 

is considered to be uncertain by a factor of 2. In L1517B, 
for example, if Ki, — 0.01 cm^ g~^, the gas will be in hy- 
drostatic equilibrium assuming Td = lOK. This higher 
value is not too far from the 0.00 8 cm'^ g^^ recommended 
by [Osscnk opf fc HenningI l(l99l for a standard gas-to- 
dust ratio of 100 and dust in regions of density w 10^ 
cm~^, which is typical of the L1517B center, and recent 
observations do in f act suggest a systematic increase of 
K„ at high densities ijKramer et al. ll20o3 iBianchi et al. I 
|2003|) . Thus, it is still possible that L1517B is exactly or 
very close to a state of hydrostatic equilibrium where its 
(mostly thermal) pressure balances its self gravity. This 
higher value, however, will not bring L1498 to a state 



^ The parameters directly derived from fitting the radial 
distribution of dust emission with an isothermal model (in- 
cluding an on-the-fly simulation) are the King radius rx = 
^9o-2/47rGpo (from the shape of the profile) and the cen- 
tral emissivity ji, = Hi^poB{Td) (from the intensity scale), 
where po is the central density and Td the dust temperature. 
Substituti ng pp from the second equation into the first one, we 
obtain k^, B{Td). 



of equilibrium, and this is also consistent with the non 
spherical shape of this core. 



4.2. Core chemistry and comparison with models 

The abundance profiles of NH3, N2H+, and the iso- 
topomers of CO and CS can be qualitatively understood 
as the result of the selective freeze out of molecules 
onto cold dust grains as the gas becomes centrally con- 
centrated during core formation. The lower binding en- 
ergy of N2 to grain surfaces makes this species rela- 
tively overabundant at high densities, and this in turn 
favors the presence of molecules composed of N and H 
(like NH3 and N2H"'") in the core interiors. This se- 
lective freeze out of molecules during core contraction 
and its resulting chemistry has been modeled in de- 
tail by a number of a u thors over the l ast few years 
llBergin fc Langer I Il997t ICharnley I Il997t ICaselli et al. I 



1999t INe iad fc Wage nblast I Il99 9': Aikawa et al. I l200ll:_ 
Li et al. i 12002; .Shematovich et al. . .2003; .Aikawa et al. I 
2003(), and their results are in general agreement with the 
observations and analysis presented here. A detailed com- 
parison between models and data, however, reveals some 
differences between observed molecular behavior and the- 
oretical expectations, and in this section we comment on 
the possible origin of these differences. 

The main disagreement between our observations and 
the published models involves the behavior of NH3 and 
N2H+ near the core centers. As seen in sections 3.3 and 
3.4, both L1498 and L1517B present the same pattern of 
constant N2H+ abundance together with a factor-of-a-few 
central NH3 increase, and this pattern seems also present 
in other starless cores (see TMCWC). Most chemical mod- 
els, however, predict that NH3 and N2H+ will behave 
in a similar fashion given that they have N2 as a com- 
mon parent. This is qualitatively correct in that the two 
species trace the same regions but the q uantitative dif- 
ferenccs noted above require explanation. I Aikawa et al. I 
{2003) have made a detailed attempt to fit observations of 
L1544 (their Fig. 8). Their best model (a delayed collapse 
with N initially molecular and a high CO binding energy) 
has an essentially constant NH3/N2H+ abundance ratio 
as a function of radius. 

It is unclear which processes cause the NH3 and N2II+ 
abundances to follow different trends at high densities and 
produce the relative enhancement of NII3 . Both have their 
origins in molecular nitrogen and hence their abundance 
ratio should be independent of the N2 abundance. Their 
formation involves reactions with IIe+ (in the case of NII3) 
and Hjj' (in the case of N2H+ ) and thus their abundance 
ratio is sensitive to the abundance ratio of these ions. 
The ionic abundances are certainly sensitive to the deple- 
tion of CO and other molecules and hence uncertainties 
in key rates determining the abundances of He+ and 
may be responsible for variations in the NII3/N2H+ abmi- 
dance ratio. Apart from this, we note that while N2II"'" is 
likely mainly destroyed by recombination with electrons 
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and reactions with CO, it is rather unclear what the main 
destruction mechanism for NH3 is in regions where CO 
(and hence carbon compounds in general) is sufficiently 
depleted that reactions with ions such as C"*" become neg- 
ligible. The answers to these questions may provide the 
explanation for the observed gradient in the NH3/N2H+ 
ratio. 

4.3. Core formation 

Although it is clear that dense cores contract from the 
lower density gas of their surrounding molecular clouds, 
how this process occurs and how lo ng it takes is st ill a 
matter of much uncerta inty (see, e.g. JShu et al. lll98^ and 
iHartmann et al. iboOll for two opposing views). This state 
of affairs results in part because it is not yet possible to 
infer a core contraction history directly from observations, 
as we lack reliable time markers to assess the evolutionary 
state of different cores. Fortunately, the time and density 
dependence of the freeze out process may offer the possibil- 
ity of shedding some light on the core contraction history 
by allowing us to correlate the molecular abundances of 
volatile species with the evolutionary state of the gas in a 
core. Here we present a first (and preliminary) application 
of these ideas to the case of the L1498 and L1517B cores. 

Both L1498 and L1517B appear highly symmetric in 
their distribution of N2H"'", NH3, and dust continuum 
(Fig. 3), which are the most reliable tracers of the dense 
core gas. L1498, especially in the density-sensitive N2H+, 
presents an evident bilateral symmetry with respect to a 
NE-SW axis, and L1517B has an almost perfect circular 
symmetry with only minor deviations. These symmetries 
in the emission of the most robust tracers most likely re- 
flect a high degree of symmetry in the underlying distribu- 
tions of matter in the cores, in agreement with their state 
close to hydrostatic equilibrium. 

In contrast with the symmetric distributions of N2H"'", 
NH3, and dust continuum, the distributions of CS and CO 
emission in L1498 and L1517B are highly asymmetric (Fig. 
3). As we have seen in sections 3.5 and 3.6, these asymme- 
tries arise from a pattern of non axisymmetric abundance 
variations that is superposed to the central abundance 
drop due to molecular freeze out. These non axisymmetric 
variations are located at radii where freeze out is already 
prevalent, so they very likely represent variations in the 
amount of freeze out, in the sense that the bright spots 
to the SE and E of L1498 and to the W of L1517B are 
places where this process is less severe and therefore the 
CO and CS abundances are higher. This interpretation is 
supported by the fact that the same regions have at the 
same time bright CO, CS, and other molecules that freeze 
out at similar densities (Tafalla et al. 2004 in preparation), 
indicating that the process responsible for the asymmetry 
is not limited to one particular chemical species. 

As molecular freeze out is very se nsitive to the time 
the g as has spent at high density fe.g.. iBer gin fc La ngerl 
the simplest explanation of the lower freeze out at 



the CO and CS bright spots is that these regions have 
somehow stayed at high density less time than other re- 
gions with similar radius (and therefore similar present 
density). (Pre-existing chemical inhomogeneities seem un- 
likely because of the small scale.) For this to happen, 
the contraction of the gas in L1498 and L1517B has to 
have been non spherically symmetric, with the gas at 
the CO and CS bright spots being "younger" (less pro- 
cessed) and therefore more recently accreted. This inter- 
pretation is consistent with the correlation between en- 
hanced CO/CS abundance and "high" velocity N2H+ and 
NH3 found in section 3.7.2 (see Fig. 13), as the younger 
portions of the core will have not yet attained perfect 
equilibrium and thus retain a fraction of their contraction 
speed. Kinematics and chemical composition, therefore, 
agree in suggesting that L1498 and L1517B have formed 
from non spherical contraction of cloud material. The time 
scale of this contraction seems to be of the order of a few 
Myr (section 3.7.2). 

The above time scale, although short, is not enough 
to discriminate between the different scenarios of star for- 
mation. On the one hand, the value is con sistent with 
a pic t ure of relatively rapid star formation ijElmegreen I 
hood iHartmann et al. I l200l|) . but on the other, it is 
also consistent with ambipolar diffusion models that start 
with an initially not highly subcritical cloud l(Li I 
ICiolek fc Basu ir200l[) . In paper 2, we will revisit the time 
scale issue based on the analysis of the freeze out process 
of a large number of molecular species. Here we note that if 
the magnetically-mediated scenario is correct, the fastest 
gas motions are expected to occur along the magnetic field 
axis, and this would imply a correlation of the observed 
abundance anisotropics with the magnetic field direction. 
Observations of the polarization direction in the submil- 
limeter continuum should be able to test this prediction. 

More conclusive than the kinematic estimate of the 
contraction time scale is the measurement of the (low) 
turbulence level inside the cores. Our observations of 
L1498 and L1517B (Fig. 10, top) appear to be inconsistent 
with typical expectations from the super sonic turbulence 
scenar io (see iMa c Low fc Klesser il 120031 for a recent re- 
view) . iBallesteros^aredes^tal. 1 1)2003^ have recently ar- 
gued that density profiles similar to those of hydrostatic 
equilibrium systems can appear fortuitously in model sim- 
ulations of supersonic turbulence, showing that the shape 
of the density profile is not a reliable indication of a core 
equilibrium state. Althoug h it is t rue that the cores in the 
iBallesteros-Paredes et al."! l(2003j) simulations can mimic 
the density profiles of systems like L1498 and L1517B, 
their internal velocity fields clearly betray their "hydro- 
static dis guise." Turbulent niqdel co res in the simulations 
of Bal lcsteros-Paredes et al.1 l(2003fl present velocity pro- 
files with internal changes larger than the speed of sound, 
and on the order of 0.5-1 km s""'^ in cases of driven tur- 
bulence (their Figs. 9-11). Our observations rule out such 
strong velocity changes by about an order of magnitude, 
and not only in the plane of the sky (from our radial pro- 
files of line center velocity, see Fig. 11), but also along any 
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given line of sight (our linewidth measurements, see Fig. 
10). 

The problem of the supersonic turbulence scenario is 
unlikely to be restricted to L1498 and L1517B, as these 
objects seem representative of the Tauru s core population 
in, fo r example, having narrow lines ('e.g.. lFuller fc Mvers I 
This problem points to the need of any model 
of core formation for producing condensations which ap- 
pear thermally dominated and extremely quiescent. How 
these close to hydrostatic equilibrium configurations can 
be achieved in the relatively short time of a few Myr seems 
to pose a challenge to most core formation models. 

5. Conclusions 

We have presented detailed models of the internal 
structure of two close-to-round starless cores in the 
nearby Taurus-Auriga star formation complex, L1498 and 
L1517B. Our models have been constrained from the si- 
multaneous fit of the spatial distribution and the central 
emerging spectrum of at least two transitions of NH3, 
N2H+, CS, C^^S, Ci«0, and C"0, together with maps 
of the 1.2 mm continuum. The main conclusions of this 
work are: 

1. Both cores present gas temperature distributions 
consistent with a constant value (10 K for L1498 and 9.5 
K for L1517B) and a very small rms (1 K). No evidence 
for radial temperature gradients is found in the central 0.1 
pc. 

2. Assuming constant dust temperature and emissivity, 
we derive density profiles for both L1498 and L1517B of 
the form n{r) = no/{l + (r/ro)"), where no, tq, and a 
are free parameters. For L1517B, the a — 2.5 best fit is 
indistinguishable from an isothermal sphere fit, and for 
L1498 the a — 3.5 best fit is close to the isothermal fit. 

3. From the combination of radiative transfer modeling 
and hyperfinc fitting of the NH3 and N2H+ spectra, wc 
determine that in both L1498 and L1517B the non thermal 
component of the linewidth is constant with radius over at 
least the central 0.1 pc. This component is smaller than 
the thermal linewidth and has a FWHM of about 0.12 
km in both cores. Its scatter in the radial profiles 
is extremely small (< 5%), and it is consistent with the 
measurement uncertainty. 

4. Using a Monte Carlo radiative transfer code and 
the radial profiles of density, temperature, and linewidth, 
we determine radial profiles of molecular abundance for 
all observed species. For NH3, the abundance at the core 
centers is a factor of a few higher than at intermediate 
radii (0.05 pc). For N2H+, a constant abundance model 
fits well the emission at all radii in L1517B, while a drop 
by a factor of 3 at large radii is preferred in L1498. For 
both CO and CS (and isotopomers), constant abundance 
models fail dramatically to fit the emission, and an abrupt 
central abundance drop is needed to explain the data. 
The abundance drop of these two species is so severe that 
the emission can be modeled with step functions and no 
molecules at the center. 



5. The central abundance drop of CO and CS seems 
to arise from the freeze out of these molecules onto cold 
dust grains at densities of a few 10'* cm~'^. The simul- 
taneous survival of N2H+ and NH3 seems related to the 
lower binding energy to grains of N2, which allows the 
production of these two species at the core centers. The 
relative increase of the NH3 abundance over N2H+ prob- 
ably results from a combination of a more efhcient NH3 
production in the highly CO-depleted core interiors and a 
decrease in the ion fraction with density, although further 
modeling is required to understand this behavior. 

6. Superposed to the pattern of radial abundance gra- 
dients there is a weaker pattern of non axially symmetric 
abundance variations in both CO and CS. This pattern 
results from asymmetries in the distribution of CO and 
CS freeze out, and its presence indicates that the process 
responsible for the freeze out has not been spherically sym- 
metric. As gas contraction to form the cores is the likely 
cause of the central freeze out, the asymmetry suggests 
that core formation has occurred in a non spherically sym- 
metric way. 

7. The analysis of the line center velocities of the freeze- 
out resistant N2H+ and NH3 reveals spatial variations at 
the level of 0.1 km s"* in both L1498 and L1517B. The 
spatial distribution of these variations is complex and can- 
not be explained with a simple pattern like rotation. It 
most likely represents internal bulk motions in the core 
gas. Some gas at anomalous velocities has a spatial dis- 
tribution similar to the gas with anomalously high CO 
and CS abundance in the pattern of non axisymmetric 
freeze out. This similarity suggests that part of the ob- 
served bulk motions result from residual core contraction. 
If so, the contraction time is of the order of 1 Myr. 

8. We have used the model physical parameters to 
study the stability of L1498 and L1517B. The small non 
thermal linewidth indicates that any turbulent pressure 
support is negligible in both cores (< 10 %), and that in 
the absence of magnetic fields, thermal pressure is the only 
stabilizing agent. Although the density profiles have the 
shape expected for isothermal equilibrium, the mass esti- 
mated from the continuum is larger than what could be 
held in equilibrium. Either our choice of dust parameters 
is incorrect, there is additional magnetic support, or the 
cores are in state of (slow) contraction. In any case, the 
observed velocity and linewidth profiles in both L1498 and 
L1517B are inconsistent with representative predictions of 
supersonic turbulence models that generate density pro- 
files similar to those of an isothermal sphere. 
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Appendix A: A simple analytic approximation to 
the isothermal sphere 

In TMCWC we fitted the volume density profiles of 5 dense 
cores using an analytical expression of the form 

n(r) — ; ; ; , 

with no, ro, and a as free parameters. The motivation of 
this choice was purely empirical, based on the ability of these 
profiles to fit the observed radial distribution of dust con- 
tinuum emission. A later hydrostatic equilibrium analysis of 
these profiles (unpublished) showed that for a = 2.5 an al- 
most perfect balance between isothermal plus constant turbu- 
lence pressure and self gravity was achieved. Intrigued by this 
"coincidence," especially because an a = 2.5 profile lacks the 
proper asymptotic behavior for an isothermal sphere (l/r'^), 
we carried out a point-by-point comparison between this pro- 
file and__th£numer^^ of the isothermal function given 
bv lChandrasekhar fc Wares I (119491) . These authors calculated 
numerically the mass density profile for an isothermal sphere 
normalized to the central density (po), which can be written 
as p — exp(— t/j) with being the solution of 



and ^ = \J iTvGpo /cr^ r (see IChandrasekhar Ill939l for a full 
discussion). To compare our a = 2.5 profile with this function, 
we normalize it and require that it reaches the half maximum 
density at the same radius as the isothermal solution (^ ~ 
2.25), and this completely fixes the profile as 

P-ff (?) = i + (^/2.25)2-5 

Fig. Al shows (top) a compariso n between the a = 2.5 pro- 
file an d the numerical solution from IChandrasekhar fc Wares I 
il949l) . For additional comparison, we al so present the classi- 
cal m odified Hubble law approximation jBinnev fc Tremaine I 
[l987) a nd the re cent simple approximation proposed by 
[Nataranian fc Lynden-BcU (1997) (their case b). As can be 
seen, the a = 2.5 profile remains within 10% from the exact so- 
luti on up to a larger (, value than the o ther solutions (although 
the lNataranian fc Lvnden-Bell I (119971) profile finally converges 
to the exact solution for large 5). In fact, the a — 2.5 solution 
is accurate within 10% for £ < 23, mor e than 3.5 times the 
instability radius dEbert lll955i r iB'onnor | [l956l. and over this 
range, the density varies by more than a factor 300. The as- 
sociated instability point of this solution ((, — 6.467) is within 
0.25% of the exact solution (^ = 6.451. lHunt'er1ll977f) . This 
shows that the a = 2.5 profile can be used to approximate any 
stable Bonnor-Ebert sphere within 10%. 

Appendix B: Kinetic temperature estimate from 
NH3 

The lack of allowed radiative transitions between the 
(J,K) = (2,2) and (1,1) levels of NH3 makes their relative pop- 
ulation (described by the rotational temperature T^^) highly 
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Fig. A.l. Top: Numerical solution of the isothermal func- 
tion from lChandrasekhar fc Wares I (| 19491) (squares) com- 
pared to our simple a = 2.5 analytic model (line). Bottom: 
relative error of several analytic approximations to the 
isothermal sphere. The solid line represents the a — 2.5 
profile presented here, the long-dashed lines is the mod- 
ified Hubble law, and the dotted line is model b from 



iNataranian fc Lvnden-Bell Note how the a = 2.5 

profile remains within 10% of the exact solution to higher 
^ values than t he other profiles (although i t fina lly di- 
verges while the lNataranian fc Lvnden-Bell I ()l997|) curve 
converges to the right solution at large ^). 



sensitive to collisions, and therefore a good estimator of the gas 
kinetic temperature. In the traditional NH3 analysis, the (1,1) 
and (2,2) populations are calculated from the spectra assuming 
LTE conditions, and the derived rotational temperature is 
converted into an estimate of t he gas kinetic temperature Tk 
(|Wahnslev fc Ungerechts Il983h . To test this procedure and to 
calibrate the T^-Tk relation using realistic radiative trans- 
fer and cloud conditions, we have run a series of Monte Carlo 
models and produced synthetic (1,1) and (2,2) spectra. These 
spectra have been analyzed using the standard NH3 procedure, 
and the derived gas parameters have been compared with those 
used in the models. 

Our NH3 Monte Carlo mo del has been described in 
TMCWC and is based on that of iBernes I (Il979t) . adapted to 
handle NH3. It includes the 6 lowest levels of para-NH3, which 
contain more than 99.9% of the molecules at 10 K (typical 
core gas temperature). It assumes that the relative hf sub- 
populations of the meta-stable levels are in LTE (as required by 
observations), and calculates the relative populations between 
the (1,1), (2,1), and (2,2) levels in a self-consistent manner. It 
uses the collision coefficients with H2 of Danbv ct al. ( 1 98^ . 
which seems to be the most accu rate set available, ac cording to 
recent laboratory measurements (IWiUev et al. Hiooi). As cloud 
parameters, we use the density profiles derived from dust emis- 
sion for L1498 (L1498-m,odels) and L1517B {L1517B-models), 
and assume a constant gas temperature (given the results of 
section 3.3), which we have varied from 5 to 20 K. 

In all 9 models we have run we find that while both the 
(1,1) and (2,2) excitation temperatures vary with radius by a 
factor of several (as a result of the core density gradient), the 

rotational temperature stays constant over the core within 
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Fig. B.l. TJ,^ versus Tk derived using Monte Carlo mod- 
els. Tk is the (constant) gas kinetic temperature and T^ 
is the (2,2)-to-(l,l) rotation temperature derived by ap- 
plying the standard NH3 analysis to the Monte Carlo spec- 
tra. Filled triangles and open squares represent models 
having the analytic density profile derived in section 3.1 
for L1517B and L1498, respectively (see Table 2). Muhiple 
points at the same Tk represent values measured at radii 
equal 0, 10, 30, 60, 100, and 150 arcseconds (data con- 
volved to the 40" resolution of the Effelsberg 100m tele- 
scope). The solid line is the analytic expression proposed 
in the text to derive Tk using the T]^ estimate from the 
standard NH3 analysis. 

2% or less. This makes the derivation of using the standard 
NH3 analysis (e.g., Bachillcr ot al. 1987) rather accurate, de- 
spite of violating one of its basic assumptions: that of a single 
excitation temperature for both (1,1) and (2,2). In fact, our 
models show that in the range of kinetic temperatures we have 
explored (5-20 K), estimates using the standard analysis 
agree within 5% with the real values, and that the agreement 
is at the 1% level for temperatures around 10 K or lower. Thus, 
T]1 can be easily and accurately determined from observations 
for realistic core conditions. 

The critical part of the NH3 analysis is the conversion of 
the estimate into an estimate of the gas kinetic temper- 
ature Tk, as this step requires a full solution of th e radiative 
transfer equation. IWalmslev fc Ungerechtsi il983l) have pre- 
sented a simple analytic expression derived under simplifying 
assumptions and using an old v ersion of the NH3-H2 collision 
coefficients from iGreen I il98ll) . In order to derive a revised 
expression using re alistic radiative trans port and the new col- 
lision coefficients of iDanbv et al. I lll98d) . we compare the Tk 
inputs in the Monte Carlo models with the T'^ values derived 
from the above analysis. As can be seen in Figure B.l., for any 
given Tk between 5 and 20K, the L1498-models and L1517B- 
models predict the same T^ within 5% or better, and this 
shows that at least for our conditions, the T^-Tk relation is 
almost independent of the details of core density, size, etc. To 
derive an analytic e xpression for this relation , we us e as inspi- 
ration equation 2 of lWalmslev fc Uneerechts I l)l983|) , although 
we look for the inverse function of theirs because we want to 
predict Tk given the T^ value determined with the standard 
NII3 analysis. After some experimentation, we have found that 



the following expression 

rp^ _ £_R 

1-^ ln[l + l.lexp(-16/rf)] 

fits the points in the range Tk = 5-20 K to better than 5%, and 
therefore provides an accurate gas temperature estimate based 
on easy-to-measure quantities. We recommend this relation to 
derive Tk in dense, cold cores at temperatures lower than or 
around 20 K (extrapolation to higher temperatures requires 
further modeling). 
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